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ABSTRACT 

We study how runaway stellar collisions in high-redshift, metal-poor star clusters 
form very massive stars (VMSs) that can directly collapse to intermediate-mass black 
holes (IMBHs). We follow the evolution of a pair of neighbouring high-redshift mini¬ 
haloes with high-resolution, cosmological hydrodynamical zoom-in simulations using 
the adaptive mesh refinement code RAMSES combined with the non-equilibrium chem¬ 
istry package KROME. The first collapsing mini-halo is assumed to enrich the central 
nuclear star cluster (NSC) of the other to a critical metallicity, sufficient for Popula¬ 
tion II (Pop. II) star formation at redshift z ss 27. Using the spatial configuration of 
the flattened, asymmetrical gas cloud forming in the core of the metal enriched halo, 
we set the initial conditions for simulations of an initially non-spherical star cluster 
with the direct summation code NBODY6 which are compared to about 2000 NBODY6 
simulations of spherical star clusters for a wide range of star cluster parameters. The 
final mass of the VMS that forms depends strongly on the initial mass and initial 
central density of the NSC. For the initial central densities suggested by our RAMSES 
simulations, VMSs with mass > 400 M 0 can form in clusters with stellar masses of 
10"^ Mq, and this can increase to well over 1000 Mq for more massive and denser 
clusters. The high probability we find for forming a VMS in these mini-haloes at such 
an early cosmic time makes collisional runaway of Pop. II star clusters a promising 
channel for producing large numbers of high-redshift IMBHs that may act as the seeds 
of supermassive black holes. 

Key words: galaxies: high-redshift ~ quasars: supermassive black holes - galaxies: 
star clusters: general. 


1 INTRODUCTION 

The most simple path to a black hole is the death of a 
massive star of mass M > 8 M©; however, these stel¬ 
lar mass black holes cannot accrete matter fast enough 
to produce the population of observed supermassive black 


the Eddington rate ( 

Volonteri et al.l 20031: 

Volonteri & Reed 

120051: iHaimanl 120061; 

Lodato & NataraianI 

20061). The nrob- 


lem would be somewhat alleviated if Population III (Pop. 
Ill) stars were hundreds of solar masses as early simula¬ 
tions predicted (|Bromm et al ][2 oo 3) . Direct collapse of these 
stars into black holes of a similar mass could then provide 
a promising route to an SMBH if the resulting black hole 
remnants could be fed continuously. Recent, higher resolu¬ 
tion simulations predict, however, that the masses of Pop. 
Ill stars are significantly lower due to fragmentation which 
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forms a small stellar association rather than an individ¬ 
ual star. This casts doubt on the possi bility that Pop. Il l 
stars produce massive black hole seeds dCreif et al.ll201ll ). 
Even if Pop. Ill stars can reach masses of « 1000 M© 
(see e.g. iHirano et al1l2014h . there are still major issues re¬ 
garding the accretion once the bl ack hole forms which are 
mainly due to radiat ive feedback l|,lohnson fc Brommll2007l : 
IPark fc R.icottilboill ). 


A different route to massive black hole seeds is the 
direct collapse of massive gaseous cores in the centres 
of atomic cooling haloes with Tvir ^ 10,000 K at high 
redshift. In primordial galaxies that remain free from 
molecular hydrogen and metals, the gas cools very inef¬ 
ficiently at lower temperatures than the atomic cooling 
limit and matter will continue to accrete on to the cen¬ 
tral object at a high rate. Simulations indicate that the 
gas may the n directly collapse into a black hole of ~ 
10"^ — 10^ llLp eb fc Rasidll994lEisenstein fc Lo^ll995l: 


u —lu i mueo <y. r\.asu.j^|iua^. [iijisensiein (V. 

hgelman et alT ^OOl : Regan fc Haehneld ^OOl : Choi et al.l 
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l2013h . This mechanism struggles, however, if the gas can ef¬ 
ficiently fragment and cool and it is essential that the halo 
remains free from metal, dust, and H 2 contamination. 

Simulations of the direct collapse scenario therefore in¬ 
voke a strong, uni form, Lyman-Werner (LW) background 
to dissociate H 2 (IWolcott-Green et al.l I 2 OI 1 I 1 . The most 
likely environment to find such a strong LW background 
is in the vicinity of a much larger host galaxy which 
has already undergone significant Pop. II or Pop. Ill star 
formation. There is some debate about the amplitude 


tion (see e.e. IShanff et al.| 

20inl: IWolcott-Green fe Haiman 

201ll:lsugimura et al.ll2014l: 

Latif et al.ll2014blld: 

Regan et al. 

2014bll especially in the presence of cosmic 

ray and X- 

ray radiation (jlnavoshi & Omukaill201ll: llnavoshi & Tanakd 


LW background is not uniform, the critical value needed to 
dissociate th e H 2 is mu c h high er than that for a uniform 
background. iLatif et al.l (l2014al l suggest that when mod¬ 
elling the radiation spectra of Pop. II stars correctly and 
including the impact of X-ray ionization, the number density 
of direct collapse black holes decreases below that required 
to grow observed high-redshift SMBHs. X-ray feedback from 
the initial gas accretion on the black hole may further limit 
how massive such an obje ct can grow in a short period of 
time llAvkutalp et al.ll2014l L Thus, it appears prudent to ex¬ 
plore other mechanisms for forming massive black hole seeds 
at high redshift. 

In this work, we study how stellar collisions in high- 
redshift, dense star clusters lead to the runaway growth 
of a single star. As stars collide, the mass and radius in¬ 
crease which boosts the probability for future collisions. This 
process becomes remarkably unstable, and analytic work, 
as well as simulations, demonstrates that, under the right 
conditions, runaway stellar collisions can produce a very 
massive star (VMS) that may collapse t o an intermediate- 
mass black hole (IMBH) of ~ 1000 ( [Beaelrnan fc Reed 
1 19781 : iPortegies Zwart et al.l l2004 iFreitag et al.ll2006l l. For 
this reason, high-redshift nuclear star clusters (NSCs) are 
very promising can d idates for the forma t ion o f IMBHs 
(lOmukai et al.luOO^ ') . iDevecchi fc Volonteril (l2009fl used an¬ 
alytical models to demonstrate how a population of IMBHs 
might form in clusters at the centres of high-redshift galax¬ 
ies. Runaway collisions in Pop. II star clusters may therefore 
be key for explaining the presence of black holes over the en¬ 
tire observed mass range. 

Here, we use a combination of hydrodynamic simula¬ 
tions and direct summation, V-body simulations to model 
collisional runaway in dense stellar clusters at high redshift. 
In Part I (Section 2), we begin with self-consistent cosmo¬ 
logical simulations performed with the RAMSES code and 
identify dense baryonic clumps in protogalaxies for which 
we can determine detailed chemo-thermodynamical prop¬ 
erties. In Part H (Section 3), we extract the clumps from 
the cosmological simulations and populate them with stars, 
varying the stellar initial mass function (IMF) as well as 
a range of other parameters, and use this as the initial 
conditions for direct A^-body simulations performed with 
NBODY 6 . Our model is nearly self-consistent, barring the 
ability to resolve the formation of individual stars within 
the cosm ological framework, which is only now becoming 
possible (ISafranek-Shrader et al.ll2014l L 


2 PART I: COSMOLOGICAL 

HYDRODYNAMIC SIMULATIONS 

2.1 Set-up of the cosmological simulations 

2.1.1 Hydrodynamic, gravity, and chemistry solver 

We use the publicly available adap tive mesh refinement 
(AMR) code RAMSES (lTevssierll2002l ~) to follow the detailed 
hydrodynamics of the first collapsing objects at high red¬ 
shift. We have replaced the default cooling module in RAM- 
SES with the non- equilibrium chemistry solver kromeQH 
dCrassi et al ][ 2 ^. KROME uses the high-order DLSODES 
solver to solve the rate equations for the chemistry network. 
We follow the detailed abundances of 12 species: H, e“, H*^, 
H", He, He+, He++, H 2 , H^, D, D+, and HD for the reac¬ 
tions fisted in the KROME react_primordiaLphotoH2 net¬ 
work. Cooling due to metals via line transitions is included 
at T < 10'‘ K for O I, C H, Si II, and Fe H. The abundances 
of these species are pinned to the hydrogen density in each 
cell assuming a metallicity which can change throughout the 
simulation. This is further described in Section 2.2.1. We em¬ 
phasize that we do not assume an amplified LW background 
that would dissociate H 2 molecules and thus prevent cooling 
below « 10“* K and inhibit gas fragmentation. 


2.1.2 Optimizing resolution and refinement 


For the purpose of our work, it is important to choose a 
large enough cosmological box to have a sufficiently massive 
halo forming within. At the same time, we need to choose 
the maximum level of refinement of the simulation such that 
we can resolve high enough densities for star formation to 
occur while at the same time making sure that no numerical 
fragmentation happens. iLada et al. llm^) suggest that star 
formation can begin to occur at vol ume number densities 
n > 10'^ cm~®. ICeverino et al.l (l2010l ') argue that the Jeans 
length, Xj, must be resolved by at least Nj = 7 cells at the 
maximum level of refinement, Imax, to prevent numerical 
fragmentation. Note that this is higher than the often used 
value of Nj = 4 cells suggested bv lTruelove et al] lll997l L To 
study the properties of the birth clouds of NSCs that form 
in our simulations, we must evolve the simulations past the 
point of first collapse which makes these simulations sus¬ 
ceptible to numerical fragmentation (iRob ertson &: Kravtsov! 
I2OO8I : ICeverino et al.l 120101 : IPrieto et al.l 2013h . In order to 
prevent this, we implement an artificial temperature floor 
based on the Jeans criteria at the maximum level of refine¬ 
ment, 


rp _ GNj pmppL 

7r7feb22*max ■ 

Here, L is the physical length of the box at the red¬ 
shift of interest, 7 is the adiabatic index, Nj is the num¬ 
ber cells we wish to resolve the Jeans length with, p is 
the mean molecular weight, p is the mass density of the 
cell, kh is the Boltzmann constant, and /max is the maxi¬ 
mum level of refinement. We set Nj = 8 , double the number 
suggested bv lTruelove et ahl lll997ll and slightly larger than 
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the value suggested bv ICeverino et al] (120101 '). The minimum 
temperature of our simulation is governed by the physical 
temperature floor set by the cosmic microwave background 
(CMB) temperature at a given redshift. Thus, as long as 
Tfloor < 2.725(1 + z) K, we will continue to accurately re¬ 
solve the chemical and hydrodynamical properties of the gas. 
Inserting Tfloor < 2.725(1-1-2) K, Nj = 8, n = lO'^ cm“®, and 
^ = 1.22 into equation (1), we see that we can safely resolve 
this density, choosing L = 500 comoving kpc h~^ at « = 30 
with /max = 19. For these values, Tfloor = 4.9 K which is 
far below the temperature. Tomb = 84.5 K, we expect the 
gas to cool to. At coarser levels < /max, we refine in order 
to resolve the Jeans length by 16 cells, double our choice of 
Nj. In addition to these criteria, we have also implemented 
refinement criteria when the number of dark matter parti¬ 
cles per cell becomes greater than 64 as well as when the 
baryons in the cell reach the equivalent scaled mass. 

We emphasize that our choice of Nj is unlikely to be 
sufficient to resolve the turbu lent properties o f the g as on 
the scales t hat w e simulate. iFederrath et al.l (l201ll i and 
iTurk et al.l ll2012f) demonstrate that in order to capture 
these properties in simulations, the Jeans length must be 
resolved by more than « 32 cells especially in the presence 
of magnetic fields. As this is not the aim of our work and 
we only look to model the general structure and mass of the 
birth cloud of a high-redshift primordial star cluster, our 
choice of resolution should be sufficient. 


2.1.3 Initial Conditions 


We use the software package MUSIC llHahn fc Abelllioill f to 
construct initial conditions for a collisionless (dark-matter- 
only) simulation using second-order Lagrangian perturba¬ 
tions on a uniform grid at 2 = 150 with 256^ parti¬ 
cles in a 500 kpc h~^ comoving box. For the cosmolog- 
ical parameters, we assum e the most recent values from 
IPlanck Collaboration et al.l (l2013f ) {h = 0.6711, flm = 
0.3175, Ua = 0.6825, as = 0.8344). The transfer function 
used to generate the init ial conditions was created using 
CAMS jLewis_et_alJ_l2000|). We use the ROCKSTAR halo 
finder (iBehroozi et al.l2013h to identify the most massive 
halo in our simulation at 2 ; = 20 which has a mass of 
Mvir = 2.48 X 10^ Mq. We define a cubic Lagrange region 
of side 127 comoving kpc to encompass all of the particles 
at 2 = 20 within a region much larger than the virial radius 
of the halo. Multiple dark-matter-only simulations were run 
until an atomic cooling halo with Tvir ^ 10'* K was identified 
at this redshift. An object of this mass is slightly overmassive 
for our choice of box size and represents an rms fluctuation 
of > 4cr indicating that it is a rare halo. Such rare haloes 
are likely to be incorporated into the most massive galaxies 
at lower redshifts and are thus likely sites to host SMBHs 
at 2 ^ 6 (ISiiacki et al.ll2009l : ICosta et al.ll2014l L 

Baryons are introduced into the initial conditions us¬ 
ing the local Lagrangian approximation at level 8 on the 
base grid, and both the dark matter and the baryons are 
initially placed at level 11 in the refinement region which 
gives an effective dark matter resolution of 2048® particles 
corresponding to particles of mass nidm = 1.08 Mq/i“*. The 
initial level of the refined region was determined so that the 
mass of the dark mater particle does not subject the refined 
cells to A'-body heating which can occur when the mass of 


the dark matter particle is much greater than the mass of the 
cells. Various initial resolutions were tested and it was found 
that level 11 provides an efficient compromise between par¬ 
ticle number and A-body heating effects which were found 
to be negligible at this level. 

Although the high-resolution dark matter particles are 
unlikely to cause spurious heating, it is possible that low- 
resolution dark matter particles may infiltrate the refine¬ 
ment region and cause heating due to their much greater 
masses. We have checked throughout our simulation for con¬ 
tamination of low-resolution dark matter particles and found 
none of these particles within the virial radius of the haloes. 
The maximum level of refinement of our simulation was de¬ 
fined in Section 2.1.2 to be /max = 19 so that the artificial 
temperature floor remains less than Tcmb at the redshifts of 
interest. This gives us a resolution of 0.95 comoving pc /i~*, 
corresponding to a physical resolution of 0.046 pc at 2 = 30. 


2 . 1.4 Identifying the birth clouds of nuclear star clusters 

We use the clump finder implemented in RAMSES to identify 
gravitationally contracting, bound clumps within the simu¬ 
lations. The clump finder is sensitive to a density threshold, 
a mass threshold, and a relevance thresholc0. To identify a 
clump, we set the minimum density to be (1 -I- 2 )® cm“® 
which corresponds to 3 x 10* cm“® at 2 = 30, and the 
relevance threshold to 1.5. We do not put a constraint on 
the minimum mass. We choose the density criteria slightly 
higher than the density we wish to resolve in order to mini¬ 
mize the chance of identifying spurious clumps at the density 
of interest. This algorithm is used only to identify the loca¬ 
tion of clumps within the simulation. In order to calculate 
clump properties, we perform a secondary analysis where we 
define the outer edge of the clump to be the radius where 
the average density drops below n = 10* cm“®. We only 
consider clumps with A > A® = 512 cells which are all 
at the highest level of refinement. For /max = 19, the mini¬ 
mum volume of a clump is then 0.05 pc® at 2 = 30 which, 
assuming a spherical structure, sets a minimum radius of 
the clump to be Adump > 0.23 pc. The Arches cluster, the 
densest known star cluster in the local Universe with a cen¬ 
tral density of « 10® Mq pc“®, has an inner core radius 
of ~ 0.2 pc which is just about resolved by our simulation 
dEspinoza et al.ll2009lj . This suggests that the full radius of 
high-density clusters will likely be sufficiently resolved out 
to their outer radii by our simulations. 

2.1.5 When to end the simulation? 

We end the simulation when high-mass stars are likely to 
form as these eventually disrupt the cluster by supernova 
feedback. We apply the following criteria to determine when 
a clump is populated with high-mass stars. For a given 
stellar IMF, dN/dM = €,{M), the total number of stars 
above a certain threshold mass, Mthresh, is given by A*(> 
Mthresh) = ^{M)dM. The main-sequence lifetime of 

stars as a function of their mass begins to flatten for stars 

® The relevance threshold is the ratio of the density peak to 
the maximum saddle density dBleuler fc Tevssieill20l3) . Clumps 
which do not break this threshold are considered noise. 
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with M > 40 Mq. For a metallicity of Z = 10““^ Zq, this 
corresponds to a main-sequence lifetime of « 5 Myr. Thus, 
we set Mthresh = 40 Mq. We can then compute the mass 
of the cluster, Mdump needed to have W(> Mthresh) > 1- 
For a Salpeter IMF (^(M) oc with Mmm = 0.1 Mq 

and Mmax = 100 Mq, we find Mdump = 1-3 x lO'^ Mq 
for Nt{> Mthresh) = 1. Once the NSC in the simulation 
reaches this fiducial mass of 1.3 x 10^ Mq, we allow the 
simulation to run for fl ag = 3.5 Myr before we extr act the 
clump, consistent with IPevecchi fc Volonteril (l2009ll . Note 
that the stellar IMF at high redshift is unknown and the 
chosen value of Mdump will change considerably based on 
the choice of IMF. However, given that our chosen value of 
Mdump is much lower than the masses of observed NSCs, 
this is a rather conservative assumption. 

2.2 Results from the cosmological simulation 

In order to model the formation of a Pop. II star cluster, 
we identify a collapsing cloud of gas in close vicinity to an¬ 
other already collapsed object. This allows for the first halo 
to undergo an episode of Pop. Ill star formation and en¬ 
rich the surrounding gas, post supernova, to a level suitable 
for forming Pop. II stars. The secondary collapsing object, 
however, must be located at a sufficient distance such that 
the radiation feedback from the first episode of Pop. Ill star 
formation does not disrupt the collapse 

2.2.1 The collapse of two mini-halos in close separation 

A halo of Mvir = 2.48 x 10^ Mq was identified at z = 20 
in the dark-matter-only simulation with comoving box size 
500 kpc h~^. The member particles were traced back to the 
initial conditions where they were centred and the simula¬ 
tion was reinitialized with baryons. Within this larger halo, 
two mini-haloes were identified with a distance such that the 
first collapsing object can enrich the second with metals. 

At the time of collapse, the first mini-halo has a virial 
mass Mvir = 2.7 x 10® Mq, a virial radius i?vir = 62 Pf0 
and a baryon fraction of 14.2 per cent. This halo begins to 
collapse at 2 = 31.6 and is the first object to collapse in the 
entire simulation volume. 

If w e assume that the ave rage mass of Pop. Ill stars is 
40 Mq dHosokawa et al.l 1201 il l, the average lifeti me of the 
syste m prior to supernova explosion is 3.9 Myr dSchaereil 
I 2 OO 2 I I. We model the chemical enrichment from these first 
supernovae by implementing a metallicity floor of 10’ 
at 2 = 30.7. While this choice of metallicity floor is arbitrary, 
a prerequisite for the formation of a Pop. II star cluster is 
that the halo be enriched to a metallicity above the critical 
metallicity for fragmentatio n. For dust-free gas, this value 
is approximately 10“^ Zq dSchneider et al.|l2012h . Slightly 
higher metallicities are unlikely to signihcantly affect the 
hydrodynamics of the gas at the densities our simulations 
probe, because as we will see, the temperature of the gas 
in the second mini-halo reaches the CMB temperature floor 
which prevents further cooling. A slightly higher metallicity 
may only accelerate this process. 

^ The virial radius is calculated assuming that the average den¬ 
sity of the object equals 200pcrit ■ 



X [kpc] 


Figure 1. Snapshot showing the surface density in a 100 pc cube 
of the two mini-halos 8.5 Myr after the collapse of the second 
object. The secondary collapsing halo is centered. 

The simulations of iRitter et al.l d2012l l model the trans¬ 
port of metals explicitly from the supernova of Pop. Ill stars 
and demonstrate that the surrounding medium can be en¬ 
riched to metallicities as high as 10~^ Zq. The metal enrich¬ 
ment extends all the way to the virial radius of the halo in 
only 8.5 Myr. The outflow is expected to be bipolar and the 
dense filaments feeding the halo should be only minimally 
enriched. The fi rst mini-halo studi ed here is less massive 
than the halo in iRitter et al.l ll2012h . We may therefore ex¬ 
pect more mixing along the filaments as the outflow should 
be less impeded by gas in the central regions of the halo in 
our simulation. 

At 2 = 28.9, the second object begins its collapse at a 
distance of 117 pc from the centre of the first mini-halo and 
slightly outside its virial radius. This occurs 12.9 Myr after 
the collapse of the first object. This leaves sufficient time 
for the first object to form Pop. Ill stars and enrich the 
surrounding material with metals. The first mini-halo has 
since grown to Mvir = 8.33 x 10® Mq with Rvir = 97 pc and 
a baryon fraction of 16 per cent. The second object is falling 
into the potential well of the first mini-halo along a dense 
filament (see Fig. 1) and the two objects will eventually 
merge. 

Assuming a population of 40 Mq Pop. Ill stars, we can 
estimate the expected intensity of the LW radiation at the 
location of the second mini-halo as, 

( 2 ) 

where J 21 is in units of 10“^^erg cm~^ s“'^ Hz“^ sr“^, Nph ~ 
QNffesc is the total number of photons emitted per second 
which escape the galaxy for a population of A^» Pop . Ill stars 
with mass of 40 Mq, Q = 2.903x 10"*® photons s“^ dSchaereil 
I 2 OO 2 II for a 40 Mq Pop. Ill star, /esc is the escape fraction, 
and hp is Planck’s constant. Assuming a maximal /esc = 1.0 
and a separation of 117 p c, we find J 21 = 9.37Af*. 

iRegan et al.l (l2014bli have demonstrated that for an 
anisotropic LW source, the critical value of J 21 needed to 
completely dissociate H 2 and keep the gas from cooling is 


logio Mq kpc’ 
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J 21 « 10 ®. At a maximal escape fraction, complete dis¬ 
sociation would require the formation of a small cluster 
of roughly 100 Pop. Ill stars with masses of the order of 
40 Mq at the centre of the first mini-halo whereas sim¬ 
ulations of these mini-haloes tend to predict small stellar 
associations of Pop. Ill stars (|Greif et al.ll201ll l. Our simu¬ 
lations do not include dust which should be present in the 
second mini-halo and will help the gas cool and catalyse the 
formation of H 2 . One-zone models predict a lower critical 
value of J 21 needed to disrupt the for mation of HD which 
is the dominant coolant below ~ 200 K llYoshida et ahllioOTl : 
IWolcott-Green fc Haimanll2^11^ . This is because the main 
formation channel of HD, H 2 -I- D+ —^ HD -(- H”*", requires 
the presence of H 2 which is lowered by the exposure to the 
meta-galactic UV background. Given the low mass of the 
first mini-halo and the initial separation of the objects when 
the secondary collapses, it is unlikely that the radiative feed¬ 
back from the first mini-halo will significantly disrupt the 
temperature and density evolution of the second. For these 
reasons, we do not include a meta-galactic UV background 
in our simulation. This approach is likely to be conserva¬ 
tive as even a mild UV backgr ound can increase the accre - 
tion rate on to a central NSC llPevecchi fc Volonteri|[^09ll . 
Denser and more massive NSCs should more efficiently un¬ 
dergo runaway stellar collisions, and by neglecting a positive 
contribution from the UV background we may underesti¬ 
mate the final masses of the VMSs. 

2.2.2 Evolution of the second mini-halo 

In Fig. 2 (left and middle panels), we plot the mass-weighted 
phase space diagram of density versus temperature of the 
second minihalo within 10 pc of the densest cell just as it 
is collapsing at 2 = 28.9 and then 5 Myr later at 2 = 28. 
In these early phases of the collapse, cooling is dominated 
by H 2 lowers the gas temperature and the mass fraction 
of HD continues to rise at higher densities. HD and metals 
significantly contribute at higher densities until either the 
CMB temperature floor or the artificial temperature floor 
inhibits further cooling. 

In the right-hand panel of Fig. 2, we show the phase 
space diagram at 14.3 Myr after the initial collapse. By this 
point in the simulation, much of the highest density gas has 
reached either Tomb or is affected by the artificial tempera¬ 
ture floor. One-zone models predict that the ratio nHD/nH 2 
peaks around the maximum density we probe in this sim¬ 
ulation before falling off steeply at higher densities. We do 
not probe these very high densities in our simulation and 
the HD abundance is forced to remain at a high value even 
though the gas should further collapse and lower the HD 
abundance. Some of the extra HD can diffuse out of the 
clump and lower the temperature of gas. This effect would 
cause an increase in fragmentation and decrease the mass of 
our central clump. With the inclusion of the CMB tempera¬ 
ture floor, this effect is somewhat mitigated because the gas 
can only cool to Tcmb. 

Because the small-scale fragmentation properties of the 
gas are subject to numerical resolution effects, we run two 
additional simulations in which we vary /max between 18 
and 20 and bracket our fiducial resolution. In Fig. 3, we plot 
the mass and density profiles of the halo computed using 
log-spaced bins centred on the centre of mass of the halo 


for three different resolutions. At the time of collapse (left¬ 
most column), the mass and density profiles are well con¬ 
verged among the three different resolutions. By 2 = 27.7, a 
secondary clump begins to emerge in the higher resolution 
runs which appears as a bump in the density profile; how¬ 
ever, the mass profile remains well converged. As the clumps 
begin to interact, we see some discrepancy between the dif¬ 
ferent resolution runs in the very central regions; however, 
in all cases, the mass profile is well converged out to 1.5 pc. 
In Appendix A, we discuss how the small-scale differences 
and the formation of secondary clumps affect our results. 

In Fig. 4, we return to the run with /max = 19 and plot 
the evolution of the mass and radius of the clump as a func¬ 
tion of time since initial collapse. Both properties increase 
linearly for « 6 Myr until a second clump appears in the 
central regions of the mini-halo. The properties of the pri¬ 
mary clump begin to oscillate as the presence of secondary 
fragments increases the average density at larger radii far¬ 
ther away from the centre. Despite this interaction, we can 
see in the top panel of Fig. 4 that the general trend is for the 
mass to increase. In order to obtain a smooth function for 
the mass of the clump as a function of time, we fit a piece- 
wise linear model to the data points extracted directly from 
the simulation and set the break in the function to be the 
point where the oscillations begin. We apply this technique 
for both the mass and the radius of the clump of interest. 
Note that the radius remains relatively constant after the 
break. 

10.8 Myr after the initial collapse, the central clump 
has grown to the threshold mass of 1.3 x 10'* Mq. The av¬ 
erage mass accretion on to the clump is Mdump ~ 6.0 x 
10~* Mq yr“*. The average mass of dark matter within the 
clump marginally decreases although this only affects the 
total mass accretion rate on to the clump by 5 per cent. 
At 14.3 Myr, which represents 0.8 Myr -f flag, the clump 
reaches a mass of 1.77 x 10* Mq within a radius of 1.25 pc. 

2.2.3 Internal clump structure 

We assume that the stars form in the highest density re¬ 
gions of the clump. We identify the densest simulations cells 
making up a certain fraction the total mass and define them 
to be star forming. For reasons described in Section 3.1.5, 
we set the star formation efficiency (SFE), e = 2/3. In Fig. 
5, we plot the spatial distribution of these cells as viewed 
along the axis with the highest magnitude of angular mo¬ 
mentum [Jz axis) for each of the three different resolution 
simulations within the star-forming radius of 1.25 pc from 
the densest cell at three different times. Tracing this region 
throughout the simulation, we see that at the time of initial 
collapse (top row of Fig. 5), all three resolutions exhibit a 
very similar spatial distribution. 5 Myr later (middle row 
of Fig. 5), the spatial distributions of the three simulations 
still agree reasonably well, but there is clear evidence that 
the higher resolution runs collapse further and we begin to 
see fragmentation. By 14.3 Myr, the distribution of the gas 
at the different refinement levels has completely diverged 
although the mass contained within the region remains rea¬ 
sonably consistent. The two higher resolution simulations 
have fragmented into multiple clumps while the /max = 18 
run contains only one object. Interactions between the dif¬ 
ferent clumps in the runs with higher level of refinement 
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Figure 2. Mass-weighted phase space diagrams of density versus temperature at the initial collapse (left), 5 Myr after the collapse 
(middle), and 14.3 Myr after the collapse (right) which is the point at which we extract the clump properties from the simulation for 
input into direct Af-body calculations. This includes all gas within 10 pc of the densest cell. The gas cools to a few hundred kelvin 
initially due to H 2 , cooling then reaches the CMB temperature floor due to HD and metal cooling. The dashed lines represent the CMB 
temperature floor and the artificial temperature floor as labelled. 



Figure 3. Enclosed mass profiles (top) and density profiles (bottom) for three different resolutions as a function of redshift. The dark 
matter, gas, and total profiles are represented by black, red, and green lines, respectively. The solid, dashed, and dotted lines represent 
Imax = 18, 19, and 20. 


have caused a significant deviation in the dynamics in the 
central region whi ch makes the structu res appear different. 
As pointed out bv lReean et al.l (l2014al l. interpreting the re¬ 
sults of AMR simulations at different refinement levels is 
non-trivial and choosing a refinement level that is either too 
high or too low can lead to misinterpretation of the overall 
dynamical evolution of collapse simulations. Our simulations 
indicate that fragmentation of the central clump is likely and 
therefore /max = 19 is the minimum resolution required in 
this study to resolve the central dynamics and this is the 
resolution we choose for further analysis. 

We hence use the spatial structure identified in the 
Imax = 19 RAMSES simulation to generate reasonably real¬ 
istic initial conditions for NBODY6 simulations of the NSCs. 


We caution here, however, that this approach is by no means 
unambiguous, and we aim to sample realistic rather than ex¬ 
act initial conditions. The resolution effects on the transition 
from the RAMSES to the NBODY6 simulations are further 
discussed in Appendix A where we demonstrate that re¬ 
gardless of which resolution is chosen between our two high- 
resolution runs, the masses of the VMSs which form remain 
consistent with each other. 

For the run with our fiducial resolution with /max = 19, 
the total volume of cells which represent the densest 2/3 of 
the total gas mass in the star-forming region is 0.109 pc®. 
In this simulation, we identify three distinct regions: a main 
central massive clump (clump 1), a less massive, off-centred 
clump (clump 2), and a further less massive, off-centred 
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Figure 4. Top. Mass of the clump in the secondary collapsing 
object as a function of time since the initial collapse. The data 
points represent the simulations and the solid lines represent the 
fitting function described in the text. The black line shows the 
dark matter mass, the red line shows the gas mass, and the green 
line shows the total mass. The horizontal, dashed, black line rep¬ 
resents the fiducial mass at which we expect the cluster to form 
high mass stars. The first vertical, dashed, black line is the time at 
which the clump reaches the fiducial mass and the second vertical, 
dashed, black line is 3.5 Myr after the first which is the point at 
which the clump mass is extracted. Bottom, Radius of the clump 
where the average density drops below 10“^ cm“® as a function 
of time since the initial collapse. The dashed, black vertical lines 
are at the same times as in the top panel and a piecewise linear 
fit to the data is also shown. 


clump (clump 3) as indicated in the bottom-middle panel of 
Fig. 5. The masses of these clumps are 9094.4 M©, 835.4 M©, 
and 170.4 M©, respectively. The structure, orientation, and 
relative velocities of these clumps are used in Part II of this 
work to generate a series of non-spherical initial conditions 
for direct N-body simulations. 


3 PART II: DIRECT N-BODY SIMULATIONS 

3.1 Setup of the NBODY6 simulations 

3.1.1 From birth cloud to star cluster 

To create initial conditions for the direct N-body simula¬ 
tions, a minimal bounding ellipsoid enclosing all of the cells 
for each of the individual regions shown in the bottom- 


middle panel of Fig. 5 is computecfl. We aim to create a set 
of star cluster initial conditions which have the same spatial 
structure and bulk velocity properties as the gas bounded 
by the three ellipsoids. The three clumps/star clusters will 
then be evolved together and allowed to interact. To con¬ 
struct flattened, ellipsoidal star clusters for each individual 
clump identified in the birth cloud, we first average the three 
primary axis lengths and then generate a spherical star clus - 
ter with this radius with MCLUSTER (iKupper et al.lf201ll l. 
For each initial sphere, the axial ratios are scaled and the 
positions of individual stars are rotated and translated ac¬ 
cording to the properties of the ellipsoids to reproduce the 
shapes and orientation of the individual clumps with respect 
to each other. The velocities of individual stars are initial¬ 
ized by choosing a value of the virial parameter Q (Q — 0.5 
represents a virialized cluster and Q < 0.5 represents dy¬ 
namically cold clusters) for the initial spherical star cluster 
prior to axis scaling and the velocities are rotated to be con¬ 
sistent with the reorientation of the positions. We do not 
attempt to scale the magnitudes of the velocities and only 
vary Q for the initial spherical cluster which determines how 
dynamically cold the initial cluster is. The bulk velocities of 
clumps 2 and 3 with respect to clump 1 are calculated di¬ 
rectly from the hydrodynamic simulation and added to the 
initial velocities of individual stars of clumps 2 and 3, but 
we do not include a net rotation for the central clump which 
may delay core collapse. The initial conditions for each of 
the individual clumps are then merged into one file to make 
a single input for NBODY6. 

The initial density profile of the spherical star clusters, 

J irior to axis sc aling, is varied between a Plummer model 
Plumme^^ll911^ and a constant-density model. We further 
use the fractal dimension option of MCLUSTER to model 
inhomogeneous systems where stars are sub-clustered to 
account for the spatial inhomogeneity in real star-forming 
systems. D — 3.0 represents a cluster with no fractaliza- 
tion and D = 1.6 represents a very clumpy distribution 
jCoodwin &: WhitworthI 120041 ') . We should note here that 
for both the Plummer and the fractal models not all stars in 
the initial conditions are placed within the bounding ellip¬ 
soid. Higher values of D effectively increase the initial vol¬ 
ume that contains the stars and therefore lower the average 
initial densities. 

The dark matter and gas not converted to stars are 
modelled as an external potential with a Plummer sphere 
and the virial radius is set to be 1.25 pc, matching that 
of the cosmological, hydrodynamic simulation. The initial 
mass of the external potential is 7524.7 M© and this mass 
increases at a rate of 6.04 x 10”“* M© yr~^ over the course of 
the simulations as measured from the RAMSES simulation 
to represent the gas and dark matter accretion on to the 
NSC. Accretion at this rather low rate affects the dynamics 
of the star cluster very little, but is, nevertheless, included 
for completeness. Note that accretion at higher rates could 
play an important role in the evolution of the cluster. 


® In addition to these three regions, there is a diffuse ring of mass 
located around the central clump and parallel to the Jz axis. The 
total mass in this ring is 5 per cent of the total mass in the main, 
central clump and we have included it in the mass we list for 
clump 1 but do not take into account its volume. 

















8 H. Katz et al. 


3.1.2 Gravity solver 

The evolution of these embedded star clusters and the for¬ 
mation of VMSs from runaway stellar collision s are modelled 
with the GPUaccelerated version of NBODY6 llAarsethll 19991 : 
iNitadori fc Aarsethll20l3l . The minimum energy conserva- 
tion requirement is set so AE/E < 10“®. Because the re¬ 
maining gas and dark matter are modelled as a background 
potential which is variable as a function of the gas accretion 
and expulsion rates, energy is not strictly conserved. 


3.1.3 Stellar evolution and metallicity effects 

Stellar evolution for in dividual stars is modelle d with the 
SSE and BSE packages dHurlev et al.l[20oO . l2002l l built into 
NBODY6. All stars begin on the zero-age main sequence 
(ZAMS). 

The lifetimes of the most massive stars are reasonably 
constant at M > 40 M©, and range from « 3.5 — 5 Myr 
with decreasing mass. These values are rather independent 
of metallicity. At this epoch in the cluster lifetime, the most 
massive stars begin to undergo supernovae, a process which 
is not modelled in our simulations. For this reason, all sim¬ 
ulations are stopped at 3.5 Myr when the first massive star 
evolves off the main sequence. Contrary to main-sequence 
lifetimes, radii of stars are heavily dependent on metallicity 
and only converge for lower mass stars. The stellar evolution 
package native to NBODY6 does not inherently sample the 
metallicity of our cl uster. _ _ 

The models of iBaraffe et al.l (|20nil l predict that stars 
with mass M « 100 Mq and metallicity Z = 10~'^Zq have 
radii of roughly half that of the lowest metallicity native to 
NBODY6. We apply a correction term to the radii to account 
for this. At M < 10 Mq, all radii are as computed for the 
lowest metallicity available in NBODY6 and for M > 10 Mq, 
a linear interpolation is applied so that 100 Mq stars have a 
radi us half that pred i cted b y the stellar evolution packages. 

IClebbeek et al.l ll2009ll demonstrated that at high 
metallicity, stellar winds become the dominant mode of 
mass-loss from a star and can significantly limit the fi¬ 
nal remnant mass of a collis ion produc t . We have imple - 
mented the mass-loss rates of IVink et al.l (Il999l . l200fll . 1200 il l 
and apply them to the collisionally produced stars with 
M > 100 Mq. The strength of the wind scales with the 
metallicity. Lower metallicity stars undergo far less mass-loss 
than their higher metallicity counterparts. For the metallic¬ 
ity considered here, stellar winds are inefficient at decreas¬ 
ing the mass of the VMS over the main-sequence lifetimes 
of these stars. However, real NSCs are likely to exhibit a 
range in metallicities and this will become important if the 
metallicity of the cluster is significantly increased. 

Another consequence of the decreased wind strength 
is the inability of the stars to unbind the gas from the 
cluster. This may decrease the number density of stars 
throughout the cluste r. Using models from STARBURST99, 
llLeitherer et al.lll999l L for a Salpeter stellar IMF with a 
maximum mass of 100 Mq and an instantaneous starburst, 
we can calculate the integrated mechanical luminosity of 
our star cluster by scaling the mass and Z®'®. Our simple 
model assumes that the mechanical luminosity is dominated 
by winds rather than radiation. By comparing this energy 
input to the binding energy of the cluster, we find that the 


energy input from the stars only becomes comparable to the 
binding energy of the cluster at « 3.5 Myr which is the life¬ 
time of the most massive stars in the cluster and the point at 
which we stop the simulations. This calculation also assumes 
that all of the mechanical luminosity couples to the gas ef¬ 
ficiently. We can, therefore, safely neglect gas expulsion in 
our simulations. 


3 . 1.4 Treatment of stellar collisions 


As the simulations begin with all stars on the ZAMS and 
are truncated after 3.5 Myr, the only collisions which can 
occur are those between two main-sequence stars. A sticky 
sphere approximation is used so that if the distance between 
the centres of the two stars is less than the sum of the radii, 
it is assumed that the stars have merged. All of the stars 
in the simulations are on the main sequence and we assume 
that when stars collide, the remnant is also a main-sequence 
star. The new star is assumed to be well mixed and t he new 
lifetime is given to be consistent with I Tout et al.l lll997ll . 
This results in a slight rejuvenation of the lifetime. If the 
collision product has M > 100 Mq, the resulting evolution 
is treated as a 100 Mq star. 

Smoothed particle hydrodynamic simulations which 
have studied the merges of main-sequence stars have shown 
that not all of the mass that enters a collision is neces¬ 
sarily retained. The mass ratio of the stars, as well as the 
orientation of the collision, ranging from head on to graz¬ 
ing, influences the amount of mass that is lost (iTrac et al.l 
I 2 OO 7 I : iDale fc DaviesI l2006l ~l . IClebbeek et al.l 1 I 2 OI 3 II demon- 
strated that the mass-loss is also slightly dependent on the 
types of stars which merge. Since all stars in our simulations 
are on the main sequence, we adopt approximated mass-loss 
rates consistent with the high-mass hal f-age main-sequence 
merger models of IClebbeek et ahl ll2013lf as follows: 


dM = min 


0.062 


M2 

0.7Mi 


,0.062 


(Ml -b M 2 ), 


(3) 


where Mi is the mass of the primary and M 2 is the mass 
of the secondary. The mass-loss is roughly constant for all 
mass ratios with M 2 /M 1 > 0.7. This is enforced in our equa¬ 
tion. We only calculate this mass-loss when the stars collide 
almost head-on such that the distance between the two cen¬ 
tres is less than half the sum of the two radii of the stars 
and when the orbital kinetic energy of the secondary star 
just prior to the collision is greater than the binding energy 


3.1.5 Star formation efficiency 

In the local Univers e the SFE, e = M*/(M* -|- Mg), is 
« 10 — 30 per cent llLada fc Ladal 1200^ ) . However, we do 
not know whether this is applicable to the hi gh-redshift 
Unive rse where the environment is very different. Dib et ahl 
1 I 2 OIIII have demonstrated that the SFE increases exponen¬ 
tially with decreasing metallicity with no relation to the 
mass of the birth cloud. The metallicity studied in their work 
is three orders of magnitude greater than the metallicity 
floor of 10“'*Zq used in this study which may suggest that 
the SFE in the NSC forming in our simulated galaxy could 
have e S> 35 per cent. Simulations of lPfalzner fc KaczmarekI 
1 I 2 OI 3 II show that, in order to reproduce characteristics of lo¬ 
cal compact clusters of similar mass to the NSC that forms 
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Figure 5. View along the axis of the densest cells which represent 2/3 of the total mass within a 1.25 pc radius of the densest cell for 
the three different resolutions (columns) and at three different times (rows). The top row shows the initial collapse, the middle row is 5 
Myr after collapse, and the bottom row is 14.3 Myr after the collapse, the point at which we extract the clump properties. Points in the 
bottom row represent the corresponding star forming region of the clump. The three clumps labelled 1, 2, and 3 in the bottom middle 
panel represent the three clumps which were extracted from the /max = 19 simulation and used to create non spherical initial conditions 
for simulation with NBODY6. The circle indicates a radius of 1.25 pc. 


in our simulation, SFEs of 60 — 70 per cent need to be as¬ 
sume^ This fi nding is further supported by the simulations 
of iFuiiil ll2014l) which demonstrate that a local SFE of more 
than 50 per cent is needed for the formation of young mas¬ 
sive clusters which have properties similar to such objects 


in the local Universe. Both observations and the simulation 
results discussed suggest that it is indeed appropriate to as¬ 
sume a rather high SFE in our star cluster simulations. 

For our fiducial model, we adopt e = 2/3 which is de¬ 
fined at the point at which we extract the clump properties 
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Figure 6. Left. Number of binaries as a function of time for a representative NBODY6 simulation with the top-heave Saltpeter IMF, with 
no initial binaries or primordial mass segregation, that results in a VMS with a mass of Mvms = 455.1 Mq after 9 separate collisions. 
The time at which the VMS undergoes a collision/merger is indicated with a vertical line where black represents a binary collision and 
red represents a hyperbolic collision. Right. The mass of the collisional runaway as a function of time for the same simulation. The solid 
vertical line is the lifetime of the most massive stars in the cluster and the point at which the simulation is stopped. The light and dark 
gray regions represent the mass ranges where PISNs and IMBHs potentially form. In each panel, the time at which the VMS underwent a 
collision/merger is indicated with a dashed vertical line where black represents a binary collision and red represents a hyperbolic collision. 


in the simulation. We should note, however, that the SFE 
is dependent on how the edge of the clump is defined. This 
is not necessarily consistent between simulations and ob¬ 
servations. We have chosen a fiducial density threshold of 
10 "^ cm“® which is higher than the densities at the edges of 
many local molecular clouds. Since observations may probe 
lower densities, the SFE which we quote will appear higher 
than the true efficiency if a fair comparison would be made 
between observations and the simulation. For example, if we 
consider a sphere of radius 10 pc around the densest cell and 
consider only gas with p > 10^ H cm“®, the effective SFE 
would be equivalent to 22 per cent. 

Because mass is accreting on to the NSC, the SFE cal¬ 
culated at the beginning of the simulation is higher than 
it would be if calculated just prior to the most massive 
stars going supernova. With the mass accretion rate of 
Afciump = 6.0 X 10“^ Mq yr“^ taken from our RAMSES 
simulation, the effective SFE at the end of the Wbody sim¬ 
ulation drops to 57 per cent, 10 per cent lower than the 
initial value. The impact of varying the SFE on our results 
is further discussed in Appendix B where we relax the as¬ 
sumption of e = 2/3. 

3.1.6 Possible outcomes of runaway collisions 

There are three possible outcomes of the cluster evolu¬ 
tion that we are interested in identifying: (1) when col¬ 
lisional runaway results in the formation of a VMS with 
M > 260 Mq, (2) when stellar collisions result in the for¬ 
mation of a pair-instability supernova (PISN) which occurs 
when mergers produce a star with 150 < M < 260 Mq, and 
(3) when collisions do not lead to efficient runaway growth 
and no star exceeds the PISN mass threshold. In order to 
sample the probability of each outcome, we generate multi¬ 
ple realizations of each set of initial conditions. 


The runaway collision process begins when the high- 
mass stars sink to the centre of the cluster and dynamically 
form binaries. Encounters with other stars perturb these bi¬ 
naries by either three- (or many-) body scattering or by bi¬ 
nary exchanges. The semimajor axis of the dominant binary 
continues to shrink as the system loses energy due to these 
encounters and if the eccentricity becomes high enough, the 
two stars merge. This process of binary capture followed 
by a merge repeats until the core evolves to sufficiently low 
density or the supernovae from the first massive stars dis¬ 
rupt the cluster. To demonstrate this process in practice, we 
show in Fig. 6 the number of binaries as a function of time 
as well as the mass evolution of the collisional runaway star 
for a representative star cluster simulation which forms a 
VMS with a mass of Mvms = 455.1 Mq after nine separate 
collisions. We indicate the times of a collision and see that 
most of these coincide with the points at which the num¬ 
ber of binaries changes. The two specific instances where 
the collision time is not related to a change in number of 
binaries occur when the VMS undergoes a hyperbolic colli¬ 
sion. The masses of the secondary objects in these types of 
collision tend to be small with the average secondary mass 
of the hyperbolic collision being 1.05 Mq compared to the 
average mass of the secondary star in binary collisions being 
56.2 Mq. Hence, the hyperbolic collisions tend to contribute 
negligible amounts of mass to the overall mass of the VMS. 


3.2 Results of the direct N-body simulations 

3.2.1 Non-spherical N-body simulations 

We begin by studying the non-spherical star clusters which 
are set up to reproduce the mass distribution obtained from 
the cosmological simulations 
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Figure 7. Initial evolution of the number density of the core as a function of time for models with D = 3.0 (left), D = 2.6 (middle), 
and D = 1.6 (right). Red lines represent models that were initialized with Plummer potentials (P) while black lines started at constant 
density (CD). Initial Q = 0.5 and Q = 0.1 are indicated by the solid and dot-dashed lines, respectively. The blue line in the left panel 
shows the evolution of a comparison spherical model (S) with D = 3.0, Q = 0.5, and /c = 0.5. 


3.2.1.1 Varying the density profile, Q and D In 

Table 1, we list the initial conditions of all of the models 
sampled in this section as well as the results. We fix the 
stellar IMF to a top-heavj0 Salpeter IMF (TH_Salp) such 
that Mmin = 1 Mq, Mmax = 100 Mq and a = —2.35. For 
each of these models, 50 realizations have been run. The 
fraction of models which produce a VMS does not exceed 
12 per cent and the most massive VMSs have masses in 
the range 275-410 Mq. The average number of collisions for 
these models is rather low. This prohibits collisional runaway 
in the majority of realizations. In Fig. 7, we plot the core 
number densitjQ for this initial set of models and see that 
it increases extremely quickly within the first 0.25 Myr and 
slowly decreases thereafter. Even models which do not form 
VMSs exhibit this generic behaviour. This number density 
rarely peaks above 10® pc“® and the final number density 
always converges to values of « 2 x 10'* pc“® at t = 3.5 Myr. 
Differences in the initial degree of fractalization have a very 
small effect on the core number density. The total num¬ 
ber of collisions increases in the models with lower values 
of D. However, this does not affect the fraction of models 
which produce VMSs. This increase is likely to be due to a 
higher initial local number density in the more fractal mod¬ 
els. There is a clear difference in the peak core number den¬ 
sity between models which are initialized with a Plummer 
profile compared to those initialized with constant density. 
The Plummer models start with higher core number densi¬ 
ties and also reach higher densities. Despite this, the models 
all converge within « 1 Myr to a core number density of 
« 10® pc“®. The models which were initialized with lower Q 
also reach higher densities than their counterparts. This is 
expected because the colder clusters are more susceptible to 
the initial collapse. The increase in number of collisions for 


the coldest clusters is very marginal (see Table 1) and once 
again, there is no increase in the fraction of VMSs that are 
produced. 


3.2.1.2 Primordial mass segregation Many observa¬ 
tions suggest that star clusters ma y form with a high de¬ 
gre e of mass segreg ation (see e.g. iGouliermis et af] 12004 
and Chen et al.ll2007 ). The simulations of Moeckel fc Clarke 
ll20lir demonstrate that during the accretion phase of gas 
on to protostars, mass segregation naturally occurs own¬ 
ing to the distribution of masses of the protostars. We do 
not simulate this phase but, if this occurs before the stars 
evolve to the ZAMS, our simulations should be initialized 
wi th some degree of mass segregation. Using the algorithm 
of iBaumgardt et al.l (l2008l) which produces clusters in virial 
equilibrium for chosen degrees of mass segregation (denoted 
A = 0 for a primordially unsegregated cluster and A = 1 for 
a fully mass segregated cluster), we introduce mass segre¬ 
gation into the initial spherical clusters prior to scaling the 
axial ratios and adjusting the positions and velocities. How¬ 
ever, the introduction of mass segregation into the models 
has a very minimal effect on the results of these simula¬ 
tions (see Table 1). The mean number of collisions remains 
roughly consistent with models without primordial mass seg¬ 
regation and there is a tendency for primordially mass seg¬ 
regated models to produce slightly more PISNs. The num¬ 
ber of VMSs that form remains the same as do the average 
masses of the VMSs. Regardless of many assumptions made 
to produce the initial conditions, the masses of the VMSs 
are robust but /vms will change with initial central density 
and this is further discussed in Appendix A where we use 
the non-spherical initial conditions taken from the /max = 20 
cosmological simulation to study the formation of collisional 
runaway in a denser cluster. 


® By “top-heavy” we mean that more mass is locked in high mass 
stars for this IMF compared to another we sample which has the 
same slope but a lower Mmin- This should not be confused with 
an IMF where most of the mass is in hi gh mass stars. 

^ The core radius is computed following ICasertano &: HutI dlQSSl l 
where the density at each particle is calculated using the five 
nearest neighbours. 


3.2.2 Spherical N-body simulations 

Our simulations taking into account the flattened and asym¬ 
metric spatial distribution of the gas in the central regions 
of the galaxy in the cosmological hydrodynamic simula¬ 
tion have succeeded in producing VMSs which may directly 
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P 

TH_Salp 

0.1 

2.6 

0.0 

2.20 

8.0 ± 4.0 

38.0 ± 8.7 

54.0 ± 10.4 

71.6 

370.5 

332.5 

3.34 

P 

TH.Salp 

0.5 

1.6 

0.0 

3.96 

8.0 ± 4.0 

38.0 ± 8.7 

54.0 ± 10.4 

82.2 

398.4 

302.6 

3.04 

P 

TH.Salp 

0.3 

1.6 

0.0 

3.34 

2.0 ± 2.0 

52.0 ± 10.2 

46.0 ± 9.6 

89.2 

312.6 

312.6 

3.14 

P 

TH.Salp 

0.1 

1.6 

0.0 

4.36 

10.0 ± 4.5 

30.0 ± 7.7 

60.0 ± 11.0 

73.5 

356.7 

316.7 

3.18 

P 

TH.Salp 

0.5 

3.0 

1.0 

1.78 

6.0 ± 3.5 

32.0 ± 8.0 

62.0 ± 11.1 

80.1 

338.8 

303.6 

3.05 

P 

TH.Salp 

0.3 

3.0 

1.0 

1.80 

4.0 ± 2.8 

48.0 ± 9.8 

48.0 ± 9.8 

82.6 

300.6 

287.5 

2.88 

P 

TH.Salp 

0.1 

3.0 

1.0 

1.74 

12.0 ± 4.9 

30.0 ± 7.7 

58.0 ± 10.8 

88.2 

339.7 

326.0 

3.27 

CD 

TH.Salp 

0.5 

3.0 

0.0 

1.56 

6.0 ± 3.5 

32.0 ± 8.0 

62.0 ± 11.1 

89.4 

302.8 

288.4 

2.89 

CD 

TH_Salp 

0.3 

3.0 

0.0 

1.98 

4.0 ± 2.8 

36.0 ± 8.5 

60.0 ± 11.0 

74.0 

275.5 

268.2 

2.69 

CD 

TH.Salp 

0.1 

3.0 

0.0 

2.06 

10.0 ± 4.5 

40.0 ± 8.9 

50.0 ± 10.0 

77.8 

411.5 

305.7 

3.07 

CD 

TH.Salp 

0.5 

2.6 

0.0 

2.22 

4.0 ± 2.8 

38.0 ± 8.7 

58.0 ± 10.8 

76.1 

320.5 

310.0 

3.11 

CD 

TH.Salp 

0.3 

2.6 

0.0 

1.56 

0.0 ± 0.0 

38.0 ± 8.7 

62.0 ± 11.1 

- 

- 

- 

- 

CD 

TH.Salp 

0.1 

2.6 

0.0 

2.04 

10.0 ± 4.5 

46.0 ± 9.6 

44.0 ± 9.4 

86.4 

332.5 

302.5 

3.04 

CD 

TH.Salp 

0.5 

1.6 

0.0 

7.90 

4.0 ± 2.8 

42.0 ± 9.2 

54.0 ± 10.4 

81.3 

310.5 

289.1 

2.90 

CD 

TH.Salp 

0.3 

1.6 

0.0 

7.12 

8.0 ± 4.0 

36.0 ± 8.5 

56.0 ± 10.6 

82.7 

428.6 

321.1 

3.22 

CD 

TH.Salp 

0.1 

1.6 

0.0 

8.22 

2.0 ± 2.0 

48.0 ± 9.8 

50.0 ± 10.0 

54.7 

278.8 

278.8 

2.80 

CD 

TH.Salp 

0.5 

3.0 

1.0 

1.80 

2.0 ± 2.0 

34.0 ± 8.2 

64.0 ± 11.3 

99.2 

358.0 

358.0 

3.59 

CD 

TH.Salp 

0.3 

3.0 

1.0 

2.20 

10.0 ± 4.5 

44.0 ± 9.4 

46.0 ± 9.6 

80.7 

355.5 

302.0 

3.03 

CD 

TH.Salp 

0.1 

3.0 

1.0 

2.00 

8.0 ± 4.0 

48.0 ± 9.8 

44.0 ± 9.4 

94.9 

320.2 

295.4 

2.96 


Table 1. pinit: the initial density profile of the cluster. “P” represents a Plummer sphere and “CD” represents a constant density model. 
IMF: the initial stellar IMF of the models. Q: virial ratio of the initial spherical star clusters. Q = 0.5 represents a spherical star cluster 
in virial equilibrium and Q < 0.5 represents a dynamically colder system. D: fractal dimension of the cluster. D = 3.0 is a smooth 
cluster and D = 1.6 represents a very fractal distribution. S: initial degree of mass segregation in the cluster. S = 0 represents a cluster 
without mass segregation and S = 1 is a fully segregated cluster. iVt-on: average number of collisions over all of the realizations. /vMSi 
/piSNi /nb: percentage of realizations which produce a VMS, PISN, or neither respectively. Error bars on these values assume Poisson 
statistics. Mgeed: average mass of the star which seeded the collisional runaway for simulations which produced a VMS. MvMS.max: mass 
of the most massive VMS which formed in all of the realizations. MyMS- average mass of the VMSs which form in the simulations which 
produced VMSs. Msmbh: the average mass of a SMBH at z = 6 assuming Eddington limited accretion with a radiative efficiency of 
e = 0 . 1 . For all models, 50 realizations are simulated. 


collapse into IMBHs. The simulated star clusters become 
spherical within a few hundred thousand years. However, 
note that this time-scale is dependent on the initial rota¬ 
tion within the sub-clumps. This has not been included in 
our simulations. Rotation tends to delay core collapse and 
should prolong the phase where the star cluster is in an 
asymmetric flattened state. While the clusters do become 
spherical rather quickly, it is important to note that the ini¬ 
tial dynamics, in particular, differ between the non-spherical 
and a corresponding spherical cluster. The evolution of the 
core number density of a spherical model which has a sim¬ 
ilar initial central density follows closely the non-spherical 
models with D — 3.0 and Q = 0.5 except for an initial 
spike which is due to the presence of a secondary clump 
in the non-spherical models (see the left-hand panel of Fig. 
7). This evolution differs from the colder and fractal non- 
spherical models in that the central number density rises 
slowly and does not decrease by 1 Myr. Regardless of this 
difference, we show in the following sections that the spe¬ 
cific evolution makes little difference for /vms or Mvms , and 
the parameters that affect the results the most are the ini¬ 
tial central density and the mass of the system. While the 


change in dynamics is interesting in its own right, it appears 
not to be of particular importance for the results of our work 
here. For that reason we explore a wider range of parameters 
for spherical star clusters and emphasize that this is likely a 
safe approximation when studying parameters such as /vms 
or Mvms- 


3.2.2.1 Setup of spherical star clu sters Observations 
of local clusters as w e ll as simula tions (iLada fc Ladal [20031 : 
iGirichidis et al.l [20111 : iBatd [201^ 1 show that star formation 
occurs deeply embedded in molecular clouds and that stars 
tend to form in the central high-density regions. To set the 
radius of our spherical star clusters for a given mass, we as¬ 
sume an isothermal density profile. We calculate the radius 
of the cluster such that the mass enclosed equals eMdump 
and set this as the radius which encloses the stars up to 
a constant factor /c. The mass within this radius is redis¬ 
tributed into a Plummer sphere and the remaining gas and 
dark matter are smoothed over a second Plummer sphere 
with a radius of the original clump multiplied by the same 
contraction factor fc. 
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For an isothermal sphere, 


M(r) = —, (4) 

-rtclump 

where Mdump = + Mg and -Rdump is the initial radius of 

the clump. We can then calculate the stellar radius as, 


R, 


.Iffclump 


ef^clump 5 


(5) 


which gives an enclosed mass profile for the stars. 


M*(r) = 


M. 




( 6 ) 


Here, M.,p = [(|f/.)" + o.,p = and U is 

an arbitrary contraction factor. 

We set the radius of the gas, Rg, equal to the origiual 
radius of the clump times fc, and smooth the mass over a 
Plummer sphere, so we find that the enclosed mass profile 
for the gas is, 


Mg{r) 


Mg,pR 

(r^+alg)!’ 


(7) 


where Mg,p = [{f|/c)^ + and ag,p = %fcRg- 

Fig. 4 shows that a small amount of dark matter is 
also present within the collapsing clump of gas at the centre 
of our mini-halo. We include this mass in the external gas 
potential of our star cluster 

In Fig. 8, we plot the average initial density profiles for 
spherical clusters with fc = 0.1, 0.2, 0.3, and 0.5 and com¬ 
pare with the spherically averaged density profiles of the 
Plummer and constant density models for our non-spherical 
star clusters discussed in Section 3.2.1. Note that the Plum¬ 
mer models used in the non-spherical simulations are sim¬ 
ilar to the models used in the spherical simulations with 
fc = 0.5. The constant-density models used in the non- 
spherical simulations are less dense in the inner regions than 
all of the spherical models considered here; however, they 
maintain a higher density at larger radii. We found, how¬ 
ever, that the difference between the Plummer and constant- 
density models does not significantly affect the probabilities 
of forming a VMS nor how massive these stars become. For 
this reason, we only investigate spherical clusters with Plum¬ 
mer density profiles. 


3.2.2.2 Varying fc, Q, and D In Table 2, we list the 
initial conditions for these spherical models where we vary 
fc, Q, and D. As we have seen in Fig. 8, fc controls the 
initial central density of the cluster and we found that for 
fixed mass, the initial central density is the most important 
parameter in determining what percentage of clusters pro¬ 
duce a VMS as well as how massive the VMS can grow (see 
Table 2). In Fig. 9, we plot the mass of the most massive 
VMS that formed as well as the fraction of clusters which 
produced a VMS as a function of initial central density. Both 
the fraction and the mass of the VMS are increasing with 
the increase in initial central density. The increase of both 
of these values is reasonably linear in log p and can be ap¬ 
proximated as 

fvMS = 42.61ogio(p*.maa;) “ 231.2 per cent, (8) 

and 

Mvms ,max ■—■ 214.1 logio(p 

*,max ) - 900.1 M©. (9) 



Figure 8. Density profiles of spherical models (dashed black) 
are compared with the spherically averaged density profiles of 
the non-spherical models. The red line indicates models with a 
Plummer density profile and the blue line is for models with a 
constant density profile and the three different lines show D =3.0, 
2.6, and 1.6. The models with D = 3.0 and 2.6 are indistinguish¬ 
able while the models with D = 1.6 have a slightly lower density 
in the inner regions of the cluster. Spherical models with decreas¬ 
ing central densities represent clusters with fc = 0.1, 0.2, 0.3, and 
0.5 as annotated, fc is the fraction by which the initial radius of 
the clump is contracted to produce these models. 



0 ^.^ 

5.0 5.5 6.0 6.5 7.0 7.5 8.0 

Logio(p.,max/MQPC~^) 


Figure 9. /vms (top) and Mvms max (bottom) as a function 
of the initial central stellar mass density for the TH_Salp IMF 
assuming 5 = 0 and 6 = 0. The data points come directly from 
the simulations and error bars are Icr Poisson errors. We over plot 
the best fit lines from equations I® & I®. 
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Spherical Models 


Initial Conditions Results 


R*,vir 

pc 

/c 

Q 

D 

Ncoll 

/vms 
per cent 

/piSN 
per cent 

/ne 

per cent 

-^seed 

M© 

Mvms, max 

Mq 

Mvms 

Mq 

MsmBH 
10® Mq 

0.071 

0.1 

0.5 

3.0 

9.08 

82.0 ± 12.8 

18.0 ±6.0 

0.0 ±0.0 

83.9 

694.5 

477.9 

4.80 

0.143 

0.2 

0.5 

3.0 

4.56 

58.0 ± 10.8 

30.0 ±7.7 

12.0 ±4.9 

86.6 

494.0 

355.6 

3.57 

0.215 

0.3 

0.5 

3.0 

2.48 

20.0 ±6.3 

55.0 ±10.5 

24.0 ±6.9 

86.3 

408.8 

336.5 

3.38 

0.358 

0.5 

0.5 

3.0 

1.12 

4.0 ±2.8 

38.0 ±8.7 

58.0 ± 10.8 

85.3 

264.1 

263.6 

2.65 

0.358 

0.5 

0.5 

2.6 

1.48 

4.0 ±2.8 

44.0 ±9.4 

52.0 ± 10.2 

92.4 

382.4 

338.2 

3.39 

0.358 

0.5 

0.5 

1.6 

3.28 

4.0 ±2.8 

50.0 ±10.0 

46.0 ±9.6 

78.5 

301.0 

297.0 

2.98 

0.358 

0.5 

0.3 

3.0 

0.92 

0.0 ±0.0 

40.0 ±8.9 

60.0 ± 11.0 

- 

- 

- 

- 

0.358 

0.5 

0.3 

2.6 

1.28 

0.0 ± 0.0 

52.0 ±10.2 

48.0 ±9.9 

- 

- 

- 

- 

0.358 

0.5 

0.3 

1.6 

3.02 

6.0 ±3.5 

38.0 ±8.7 

56.0 ± 10.6 

84.6 

303.1 

326.3 

3.27 


Table 2. All models assume a TH_Salp IMF with Mmin = 1 Mq , Mmax = 100 Mq , e = 2 / 3 , and M » = 1.01 X 10 ^ Mq . virial 

radius of the stellar component, fc'. contraction factor. For all models, 50 realizations are used. See Table[T]for definitions of other values. 


Initial Mass Functions 


Name 

Mcnin 

Mq 

Mmax 

Mq 

a 

Salp 

0.1 

100.0 

-2.35 

TH.Salp 

1.0 

100.0 

-2.35 

Kroupa 

0.1 

0.5 

0.5 

100.0 

-1.3 

-2.30 

Flat 

1.0 

100.0 

-2.00 


Table 3. Table of the different IMFs sampled for the simulations 
of spherical star clusters. 

Decreasing Q from 0.5 to 0.3 does not affect /vms or 
iVcoii, similarly to the non-spherical models. We found a 
small increase in iVcoii for the more inhomogeneous mod¬ 
els similar to what we found for our non-spherical NBODY6 
simulations. This is likely to be due to the higher local initial 
densities. The value of D, however, does not change /vms- 
Based on the results in Table 2, a spherical cluster with fc 
between 0.5 and 0.3, Q — 0.5 and D = 3.0 should give very 
similar results to the non-spherical models. 

3.2.2.3 Varying the IMF As the stellar IMF is highly 
uncertain, especially at high redshift, we test the effect of 
varying the IMF. The parameters of the different IMFs im¬ 
plemented are listed in Table 3. In Table 4, we list the initial 
conditions and results of the models where the IMF is var¬ 
ied. For these models, we have assumed fc = 0.2 which was 
previously determined to have a high percentage of models 
which produce a VMS (see Table 2). 

As the mean mass of the IMF is increased, the frac¬ 
tion of realizations which produce a VMS remain consistent 
within error bars and therefore, regardless of the IMF, the 
likelihood of producing a VMS in a high-redshift NSC is the 
same for an NSC of this mass. The average number of stel¬ 
lar collisions increases, however, for models with lower av¬ 
erage stellar mass (m) due to the higher number density of 
stars. Furthermore, MvMS.max and Mseed tend to increase as 
the IMF becomes more top heavy. This is expected because 
there are simply more high-mass stars in the runs with more 
top-heavy IMFs and therefore the mean mass per collision 


will increase. We can also see that Mseed slightly increases 
with increasing m which also reflects the availability of more 
high-mass stars. 


3.2.2.4 Primordial binaries and mass segregation: 

Observations suggest the presence of a significant fraction 
of bin aries in young stellar clusters, esp ecially for high-mass 
stars jHut et al.lll992l : rSana et al.ll2009l) . A large binary frac¬ 
tion at the centres of star clusters can effectively heat the 
cluster and prevent core collapse. This decreases the max¬ 
imum value of the central density which may inhibit the 
formation of a VMS due to collisional runaway. When clus¬ 
ters are initialized with binaries, massive stars with mass 
greater than 5 Mq are preferentially put in binaries. For 
low- mass binaries, we use a period distribution consistent 
with iKroupal lll995ll and for binary stars which have a pri¬ 
mary star with mass greater than 5 Mq , we a dopt the period 
distribution consistent with ISana fc Evanj (1201 il l. We test 
the effects of including primordial binaries as well as primor¬ 
dial mass segregation for these spherical clusters and list the 
initial conditions in Table 5. 

Introducing mass segregation for models with TH_Salp 
and Kroupa IMFs makes no difference to Mvms- The same 
is true when binaries are introduced. In all cases, the av¬ 
erage VMS undergoes ~ 1 — 3 collisions after the PISN 
mass threshold and with these low number statistics, it is 
unsurprising that Mvms is roughly the same regardless of 
variations to these initial parameters. We do, however, see 
a drastic increase in Acoii when the clusters are initialized 
with large binary fractions. The increase of the total number 
of collisions in the cluster is clearly not affecting Mvms • For 
stars that are in binaries from the beginning, interactions 
with other stars can efficiently remove energy from the bi¬ 
nary and eventually drive the binary to a merger. Without 
binaries, the massive stars first have to sink to the centre 
of the cluster and dynamically form binaries before merging 
can take place. The average mass of the mergers is, how¬ 
ever, much lower for simulations with primordial binaries 
and these newly formed, merged stars tend to have very lit¬ 
tle impact on the formation of a VMS. Despite the evolution 
of the star cluster being significantly different for simulations 
with primordial mass segregation and binaries, we can make 
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Spherical Models: Varying the Stellar IMF 


Initial Conditions Results 


IMF 

-^max 

M© 

m 

Mq 

Ncoll 

/vMS 
per cent 

/PISN 
per cent 

/ne 

per cent 

-^seed 

M© 

Mvms, max 
Mq 

Mvms 

Mq 

Msmbh 
10 ® Mq 

Salp 

79.5 

0.35 

8.95 

40.0 ± 14.1 

50.0 ± 15.8 

10.0 ± 7.1 

75.2 

353.1 

295.3 

2.96 

TH.Salp 

90.4 

3.09 

4.56 

58.0 ± 10.8 

30.0 ± 7.7 

12.0 ± 4.9 

86.6 

494.0 

355.6 

3.57 

Kroupa 

88.8 

0.64 

5.80 

44.0 ± 9.4 

38.0 ± 8.7 

18.0 ± 6.0 

84.0 

433.2 

333.9 

3.35 

Flat 

96.1 

4.65 

3.44 

48.0 ± 9.8 

40.0 ± 8.9 

12.0 ± 4.9 

89.9 

591.5 

376.7 

3.78 


Table 4. Mmax: the average maximum mass star over all realizations of the cluster, fh: the average stellar mass over all realizations of 
the cluster. See Table [T] for definitions of other values. 50 realizations were run for all models except the model with the Salp IMF where 
only 20 simulations are used. 


Spherical Models: Primordial Binaries and Mass Segregation 


Initial Conditions Results 


IMF 

A^max 

M© 

m 

Mq 

S 

6 

-^coll 

/vMS 
per cent 

/piSN 
per cent 

/ne 

per cent 

-^seed 

M© 

-^VMS,max 

M© 

Mvms 

Mq 

Msmbh 
10 ® Mq 

TH.Salp 

90.4 

3.09 

0.0 

0.0 

4.56 

58.0 ± 10.8 

30.0 ± 7.7 

12.0 ± 4.9 

86.6 

494.0 

355.6 

3.57 

TH.Salp 

90.4 

3.11 

0.5 

0.0 

4.78 

66.0 ± 11.4 

26.0 ± 7.2 

8.0 ± 4.0 

85.8 

516.8 

352.7 

3.54 

TH.Salp 

90.6 

3.08 

1.0 

0.0 

3.64 

36.0 ± 8.5 

48.0 ± 9.8 

16.0 ± 5.7 

83.2 

406.1 

318.5 

3.20 

TH.Salp 

92.6 

3.10 

0.0 

0.5 

17.6 

32.0 ± 8.0 

50.0 ± 10.0 

18.0 ± 6.0 

81.9 

438.4 

319.0 

3.20 

TH.Salp 

91.5 

3.09 

0.0 

1.0 

22.0 

22.0 ± 6.6 

50.0 ± 10.0 

28.0 ± 7.5 

79.0 

445.0 

327.0 

3.28 

TH.Salp 

91.4 

3.09 

0.5 

0.5 

18.0 

42.0 ± 9.2 

44.0 ± 9.4 

14.0 ± 5.3 

81.5 

515.1 

344.0 

3.45 

TH_Salp 

90.6 

3.09 

1.0 

0.5 

16.8 

24.0 ± 6.9 

36.0 ± 8.5 

40.0 ± 8.9 

83.2 

477.7 

336.7 

3.38 

TH.Salp 

90.0 

3.09 

0.5 

1.0 

23.7 

46.0 ± 9.6 

36.0 ± 8.5 

18.0 ± 6.0 

80.1 

465.3 

353.7 

3.55 

TH.Salp 

91.2 

3.08 

1.0 

1.0 

22.7 

26.0 ± 7.2 

38.0 ± 8.7 

26.0 ± 7.2 

82.2 

595.6 

365.9 

3.67 

Kroupa 

88.8 

0.64 

0.0 

0.0 

5.80 

44.0 ± 9.4 

38.0 ± 8.7 

18.0 ± 6.0 

84.0 

433.2 

333.9 

3.35 

Kroupa 

87.2 

0.64 

0.5 

0.0 

6.18 

42.0 ± 9.2 

50.0 ± 10.0 

8.0 ± 4.0 

80.0 

410.7 

336.6 

3.38 

Kroupa 

85.8 

0.64 

1.0 

0.0 

5.24 

52.0 ± 10.2 

38.0 ± 8.7 

10.0 ± 4.5 

78.9 

435.5 

327.7 

3.29 


Table 5. Note that the TH_Salp and Kroupa models with S = 0 and 6 = 0 are the same models listed in Table |4] and are shown here 
for ease of comparison. See Table[T]for definitions of values. For all models, 50 realizations are run. 


robust predictions for the average mass of a VMS that will 
form. 

Although Mvms is only weakly dependent on most as¬ 
sumptions we have made, we do find a factor of 3 change in 
the fraction of models which produce a VMS when including 
mass segregation and binaries. While this does not affect the 
mass function of these objects, it does affect their number 
density. Introducing binaries tends to decrease the number 
of VMSs which form while no trend is evident for increas¬ 
ing the primordial mass segregation. Including the binaries 
tends to heat the core and eject the low-mass stars, and this 
influences the evolution of the core of the cluster where the 
collisions occur. The effect of changing these assumptions is, 
however, much less significant than changes to the central 
mass density of the star cluster. 

3.2.2.5 Varying the mass of the NSC Thus far, our 
work has focused on one simulation of one mini-halo with 
one NSC with a particular mass motivated by our cosmolog¬ 
ical zoom-in simulation. The Universe likely exhibits a range 
of NSC masses with varying initial central densities and we 
also want to explore how our results depend on the assumed 
mass of the NSC. Although the first collapsing halo in our 
simulation is not be metal enriched, to get an idea of the 
variance in the mass of the NSCs in the early Universe, we 
can apply the same criteria that we used to define the NSC 


in the secondary collapsing halo to determine an NSC mass 
for the first halo. The total mass of the NSC in this halo 
including gas and dark matter is « 2.3 x 10“* Mq, which 
is clearly more massive than the NSC in the secondary col¬ 
lapsing object. We run a few additional simulations where 
we multiply the total mass of the spherical fc = 0.2 model 
by some factor fm and correspondingly increase the radius 
of this model by fRR' in order to determine how massive 
an IMBH may become for different NSC masses. In Table 
6 , we list the initial conditions and results of these models 
which have been averaged over 25 realizations. For all mod¬ 
els, we assume a TH_Salp IMF and Plummer model initial 
condition^ 

It is clear that for these more massive NSCs, the major¬ 
ity are likely to form an IMBH with some producing IMBHs 
with masses greater than 1000 Mq. The average masses of 
the VMS are also significantly higher. 

3.3 Caveats and Limitations 

The initial conditions of our models, in particular the central 
number and mass densities, play a major role in determin¬ 
ing the fraction of clusters which produce a VMS. Owning 

® Note that models with larger maintain a higher average 
density at larger radii. 
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Spherical Models: Varying the Mass of the NSC 


Initial Conditions Results 


fm 

-^^max 

Mq 

m 

Mq 

S 

b 

Ycoll 

/vMS 
per cent 

/PISN 
per cent 

/ne 

per cent 

^^seed 

Mq 

VfvMS.max 
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860.9 
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Table 6. fm- fraction by which the mass of the NSC is scaled. The radius is scaled by /, 


1/3 


See Table [J for definitions of other values. 


to numerical limitations, to the best of our knowledge, there 
is no study to date which has produced an embedded cluster 
such as the one presented in this paper that had an IMF con¬ 
sistent with observations which has been evolved from the 
protostellar accretion phase all the way to the onset of the 
ZAMS. There is thus a clear disconnect between the types of 
simulations which can form protostellar cores directly from 
molecular clouds and the direct N-body simulations which 
can follow the subsequent evolution of these stars once the 
cluster has completely formed. Simulations which follow the 
initial formation and gas accretion on to protostellar cores 
explore how the physical properties of the collapsing molec¬ 
ular cloud and the accreti on dynamics affect the initial con- 
ditions of t he star cluster (l Batelll998l: iBate fc Bonnellll2005l : 
lBatell2012l : lKrumholz et'^l2012l ). We are unable to capture 
this direct link between the properties of the birth cloud and 
the initial conditions of the star cluster with our direct N- 
body calculations. In order to properly model this system, 
we would have to self-consistently predict the distribution 
function and IMF of stars from the cosmological simula¬ 
tion and run the simulation far beyond the initial stages of 
collapse. This is, unfortunately, beyond current numerical 
capabilities. We have, however, mitigated this by studying a 
large number of possible initial conditions for the star clus¬ 
ters we have simulated. We believe that we can therefore 
confidently conclude that for our choice of physically moti¬ 
vated and reasonable assumptions, our simulations predict 
collisional runaway to be a promising mechanism for the 
formation of VMSs. 


We have allowed the stars to grow far beyond the masses 
at which stellar evolution is well understood due to the lack 
of observational constraints. In this regime, it is not ob¬ 
vious whether such an object is stable and evolves like a 
normal star especially because it is created by merging. For 
low-metallicity stars, the main-sequence lifetime and radius, 
which are the most important stellar parameters for this 
work, begin to flatten at high masses. For these reasons, 
we have chosen to evolve them as 100 Mq stars. Two other 
processes along these lines are neglected in our simulations 
that will affect stellar collisions, the inflation of the col¬ 
lisional remnant’s radius and the possible increase in the 
remnant’s main-sequence lifetime. When stars undergo col¬ 
lisions, some of the kinetic energy from the collision is ab¬ 
sorbed into the envelope of the primary which can cause 
the merger remnant to have a much greater radius than 
a comparable star with that mass which ev olves normally 
on the main sequence dPale fc Davies! [ 2 OO 6 I I. The lifetime 
of this stage is likely much shorter than the main-sequence 
lifetime of the star. However, the gravitationally focusing 


cross-section of the star scales with radius so the probabil¬ 
ity of u ndergoing a collision c an greatly increase during this 
period dPale fc Daviedbood l. Further collisions will cause 
a similar effect, and thus, if the time between collisions is 
short, the inflated radius can be sustained for long periods 
of time. The mixing of the two colliding stars can introduce 
a fresh source of hydrogen into the core and can increase the 
main-sequence lifetime of the remnant compared to a star 
with a similar mass that evo l ves n ormally on the main se¬ 
quence dGlebbeek et al.lBoOSl . 1201^ 1. The prolonged lifetime 
will increase the number of stellar collisions that remnant 
might undergo. Both effects tend to improve the prospect of 
forming a VMS and by neglecting them, our results should 
be conservative in this regard. 


4 IMPLICATIONS FOR THE (EARLY) 

GROWTH OF SUPERMASSIVE BLACK 

HOLES 

We have demonstrated that high-redshift, dense, metal-poor 
NSCs are likely to host runaway stellar collisions which pro¬ 
duce VMSs and that this process is robust to a wide variety 
of assumptions. We have, however, not addressed here what 
happens to the VMS once it forms and there are also no 
observ ational constraints for stars of this mass. iHeeer et al.l 
(l2003tl predict that VMSs with the mass and metallicity as 
in our work end their lives by directly collapsing into an 
IMBH with minimal mass-loss. However, these predictions 
need probably to be considered with a healthy scepticism 
given the lack of observational constraints. Furthermore, we 
have assumed that when the stars merge, the remnant be¬ 
comes an ordinary main-sequence star. This is also uncer¬ 
tain as the impact of the merge certainly disrupts the star at 
least temporarily and the subsequent stellar evolution also 
remains uncertain. With these caveats in mind, we now as¬ 
sess briefly whether this mechanism can be responsible for 
producing the population of bright quasars at 2 « 6 — 7 
which has an observe d lower limit of 1.1 x 10“® Mpc“® 
dVenemans et al.]l2013l l. To do this, we estimate the mass 
to which IMBHs can grow by these lower redshifts and their 
expected space density. 

Under the assumption of Eddington-limited accre¬ 
tion, the mass of a black hole increases as M = 
wh ere e is the radiat ive efficiency and 
fEdd = crTc/(47rGmp) (iFrank et al.ll2003 l. The size of £ is 
somewhat uncertain and depends on the properties of the 
accretion disc which surrounds the black hole as well as the 
spin of the black hole. Our direct N-body simulations are run 
for 3.5 Myr after the point at which we extract the clump 
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from the simulation. This corresponds to the lifetimes of the 
most massive stars. By this point, the hydrodynamic simula¬ 
tion would have evolved to z = 25.88. Assuming a canonical 
value of e = 0.1 we calculate the expected mass of the black 
hole at z = 6 and tabulate the results in Tables 1, 2, 4, 
5, and 6. It is clear that regardless of many of the initial 
assumptions of the IMF, binary fraction, initial mass segre¬ 
gation, SFE, and density profile, the average masses of the 
VMSs that do form are sufficient to grow to the masses of 
the black holes powering high-redshift quasars that we ob¬ 
serve at 2 > 6 if the black holes can be continuously fed at 
the Eddington accretion rate. 

The mass range that we predict for our VMSs is well 
within th e range predicted by some simulations of Pop. 
Ill stars jHirano et al.l[2014l l. For simple models assuming 
Eddington-limited accretion, there is thus little difference 
between assuming that the growth of SMBHs starts with 
the remnants of Pop. Ill stars and assuming that the growth 
is seeded by the VMS resulting from runaway growth in 
Pop. II star cluster. What makes the latter route perhaps 
more promising is, however, the very different environment 
in which the collision al runaway VMSs form co mpared to 
that of Pop. Ill stars. I.Iohnson fc BrommI (|2007l 'l predict a 
long time delay for efficient accretion on to a black hole in 
the mass range studied in this work due to radiative feed¬ 
back on the surrounding medium. Only in haloes where the 
gas density remains sufficiently high can this be avoided al¬ 
though even if the gas remains at high densities, the radia¬ 
tive fe edback still may limit the accretion to sub-Eddington 
levels dMilosavlievic et al.ll2009l l. Mergers with surrounding 
mini-haloes may provide the dense gas supply necessary to 
overcome most of the radiative feedback and efficiently feed 
the black hole. Recall that the formation of a Pop. II NSC 
at these very high redshifts we study requires a close neigh¬ 
bouring halo which can pollute the halo in which the IMBH 
forms with metals. The haloes in our hydrodynamic simula¬ 
tions merge shortly after the IMBH is likely to form. Thus, 
it appears unlikely that an IMBH which may form in our 
simulations would suffer a long period of inefficient accre¬ 
tion. 

A further difference between a Pop. Ill remnant black 
hole and the IMBHs forming from runaway stellar collisions 
is the pr esence of a surrounding star c luster in the latter 
scenario. [Alexander Sz NataraiarJ (|2014| ~) discuss the evolu¬ 
tion of low-mass black holes with ~ 10 M© embedded in a 
star cluster fed by dense gas flows and accreting at super- 
Eddington rates due to random motions within the clus¬ 
ter. If high-redshift, dense star clusters have very top-heavy 
IMFs, this may indeed be relevant for more massive black 
hole seeds. A simil ar mechanism was indeed suggested by 
[Davies et al.l ([201lj ~l where dense inflows of gas can initiate 
core collapse and cause efficient merging in the central re¬ 
gions of the cluster. The presence of the cluster after the seed 
has formed may be key to allowing it to grow efficiently. 

Assuming that these SMBH seeds can accrete at the 
required rates, we can calculate an approximate upper limit 
on the number density of SMBHs by calculating the total 
number density of haloes enriched above the critical metal- 
licity at the final redshift, 2 /i„, such that IMBHs can grow 
to 10® M© by 2 = 6. Assuming that the average mass of an 
IMBH that forms in our simulation is 300 M©, Zfin = 20. 
The number density of SMBHs at 2 = 6 then is approxi¬ 


mately 


nSMBH — fpfffedd'^gad^ Afthresh); 


( 10 ) 


where fp is the fraction of haloes at 2/in that have been 
metal enriched above Zcrit, // is the fraction of these haloes 
that form an IMBH, /edd is the fraction of the black holes 
which can sustain Eddington-limited accretion, and ngai{> 
A/thresh) is the total number density of galaxies above the 
mass threshold, Mthresh, which are possible sites for form¬ 
ing an IMBH. We set Afthresh = 5 x 10® M©, and using 
the Jenkins mass function ([Jenkins et al.[[200l[ l. we find, 
ngai ~ 8 Mpc“®. For the mass of the NSC which formed in 
our simulation, // is likely < 0.1. There is little constraint 
on the value of fp as we have little knowledge of the metal 
enrichment of the intergalactic medium at these extremely 
high redshifts. Metal enrichment in the early Universe is 
almost certainly very patchy and confined to overdense re¬ 
gions. However, even for // = 0.01, a rather small value of 
fp « 10“® would be sufhcient to explain the observed SMBH 
number density based on this simple estimate with /edd = 1- 
The value of /edd also remains highly uncertain and it is un¬ 
likely to be unity given the plethora of environmental feed¬ 
back mechanisms that inhibit eff i cient accretion on to black 
holes ( Johnson fc BromrrJ [2007[ : [Milosavlievic et al.l [2009[ : 
[Park fc Ricottjboilh . For this to occur, black holes would 
require a constant supply of cold, low-angular-momentum 
gas, and since we do not follow the hydrodynamical simu¬ 
lation past the formation of the NSC birth cloud, it is un¬ 
certain whether such a reservoir is available. Note, however, 
that the mass function of IMBHs forming in high-redshift 
NSCs is unlikely to be a delta function at the minimum 
mass required to form an SMBH by 2 = 6 so /edd need not 
be 1 for this process to produce the population of observed 
SMBHs at 2 = 6. In addition, other mechanisms may al¬ 
low black holes to accrete at super-Ed dington rates when 
a surrounding star cluster i s present ([Davies et al.[ [201 1[ : 
[Alexander fc NataraiarJ[2014[ l. 


5 CONCLUSIONS 

Using a combination of high-resolution, hydrodynamic cos¬ 
mological zoom-in simulations and spherical and non- 
spherical direct N-body simulations, we have demonstrated 
that stellar runaway growth at the centre of nuclear Pop. 
II star clusters in high-redshift protogalaxies (20 ^ z ^ 30) 
metal enriched by nearby companions is a promising route to 
form VMSs with masses as high as 300 — 1000 M©. We find 
that the average masses of the VMSs that are produced in 
NSCs with an expected total mass of typically Af* « 10"* M© 
are relatively robust to changes in the stellar IMF, number 
of primordial binaries, initial degree of mass segregation, as 
well as initial density profile, but increase strongly with the 
increase in initial central density and total mass of the star 
cluster. If the VMSs formed in this way can directly col¬ 
lapse to IMBHs with moderate mass-loss, they are promis¬ 
ing seeds for growth into the billion solar-mass black holes 
observed at 2 « 6. Our simulations further predict an en¬ 
hanced number of PISNs as the simulations which fail to 
produce a VMS often host at least one or two high-mass 
collisions. This may result in a rapid early enrichment of 
the IGM with metals and may cause the early pollution of 
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other protogalaxies making them also susceptible to the col- 
lisional runaway process. Modelling the evolution of the host 
galaxy post-supernova will be needed to predict the effects 
of these PISNs on their environment 

The presence of a large numbers of accreting IMBHs 
as well as an enhanced rate of PISN should signihcantly 
alter the early evolution of galaxies and will have impor¬ 
tant implications for the interpretation of observations of 
the high-redshift Universe. 
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APPENDIX A: CONVERGENCE TESTS FOR 
NON-SPHERICAL INITIAL CONDITIONS 

In Fig. 3, we have shown that despite the fragmentation 
present in higher resolution simulations, the mass contained 
in the central region of the central collapsing galaxy remains 
well converged. To improve on this point, we plot the phase 
space diagrams of density versus temperature at three dif¬ 
ferent times for the Imax = 20 simulation in Fig. Al. The 
evolution is nearly identical to the Imax = 19 run which is 
shown in Fig. 2, and the main difference between these two 
simulations is caused by clump-clump interactions which af¬ 
fect the central structure. The structure of the cells we iden¬ 
tify as star forming in the /max = 20 simulation changes and 
becomes denser compared to the /max = 19 run (see Fig. 5). 
We can test how these denser clumps affect the formation 
of VMSs by applying the same method for creating initial 
conditions as used in Section 3.1.1 to the highest resolution 
simulation. 

In Table Al, we list the initial conditions and results of 
direct V-body simulation run using as input the /max = 20 
run. Even though these runs have slightly less mass, we hnd 
that /vMS has increased due to the increase in density. Over¬ 
all, Mvms remains reasonably consistent with our previous 
results which suggest that this mass is robust. 


APPENDIX B: VARYING THE STAR 
FORMATION EFFICIENCY 

Most of our models were based on the assumption of an 
SFE of e = 2/3 which is motivated by a series of obser¬ 
vations and simulations of relevant environments. We test 
two models where the SFE is changed and these models are 
listed in Table Bl. The cluster parameters are derived using 
equations (4)-(7) and note that the initial central densities 
of the clusters are dependent on the chosen e. We see that 
VMSs can form in clusters with lower M* although they 
likely require higher central densities. 
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Figure Al. Mass-weighted phase space diagrams of density versus temperature at the initial collapse (left), 5 Myr after the collapse 
(middle), and 14.3 Myr after collapse (right). This includes all gas within 10 pc of the densest cell. The gas cools to a few hundred kelvin 
initially due to H 2 cooling then reaches the CMB temperature floor due to HD and metal cooling. The dashed lines represent the CMB 
temperature floor and the artificial temperature floor as labelled. 


Non-Spherical Models: Convergence 


Initial Conditions Resnlts 


Pinit 

IMF 

Q 

D 

S 

Vcoll 

/vMS 

/PISN 

/ne 

-^seed 

M0 

MvMS.max 

Mo 

MvmS 

Mo 

MsmBH 
10® Mo 

CD 

TH.Salp 

0.5 

3.0 

0.0 

2.38 

16.0 ±5.7 

46.0 ±9.6 

38.0 ±8.7 

87.6 

311.8 

290.4 

2.91 

CD 

TH.Salp 

0.3 

3.0 

0.0 

2.40 

18.0 ±6.0 

56.0 ± 10.6 

26.0 ±7.2 

84.8 

351.3 

294.3 

2.95 

CD 

TH.Salp 

0.1 

3.0 

0.0 

2.44 

24.0 ±6.9 

32.0 ±8.0 

44.0 ±9.4 

89.0 

349.0 

296.6 

2.98 


Table Al. 50 Realizations of each model are generated. See Table[T]for definitions of the listed quantities. 


APPENDIX C: ALTERNATIVE MODEL FOR 
MASS SEGREGATION 

We generate one additional model in order to compare an 
alternative method for primordial mass segregation with the 
one we have used in the text in order to determine how 
sensitive our res ults are to the algo rithm used. Here we adopt 
the algorithm of ISubr et al.l ll2008l ') and the results are listed 
in Table Cl. For this model, we use the TH_Salp IMF and it 
is directly comparable with the model in Section 3.2 with the 
same IMF and initial degree of primordial mass segregation. 
Within the la errors on the percentages of VMSs and PISNs 
that form, these models are entirely consistent suggesting 
that our choice of algorithm is not significantly affecting the 
results. 


Log,o(M/M,o,) 
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Spherical Models: Varying the Star Formation Efficiency 


Initial Conditions Results 


e 

© 

o 

R*,vir 

pc 

/c 

Q 

D 

Ncoll 

/vMS 
per cent 

/PISN 
per cent 

/ne 

per cent 

-^seed 

M0 

A7vMS,max 

Mo 

MvmS 

Mo 

MsmBH 
10® Mq 

1/2 

0.76 

0.107 

0.2 

0.5 

3.0 

3.90 

50.0 ±10.0 

40.0 ±8.9 

10.0 ±4.5 

84.1 

492.8 

345.8 

3.47 

1 

1.52 

0.215 

0.2 

0.5 

3.0 

4.40 

45.0 ±9.5 

40.0 ±8.9 

15.0 ±5.5 

82.8 

551.0 

402.5 

4.04 


Table Bl. 20 Realizations of each model are generated. See Table^for definitions of the listed quantities. 


Spherical Models: Alternative Mass Segregation 


Initial Conditions Results 


Profile 

Parameters 

.^coll 

/VMS 
per cent 

/piSN 
per cent 

/ne 

per cent 

^^seed 

M© 

MvMS.max 

Mo 

MvmS 

Mo 

MsmBH 
10® Mo 

Subr 

S = 0.5 

3.92 

50.0 ± 10.0 

36. ±8.50 

14.0 ±5.3 

84.4 

559.5 

350.3 

3.51 


Table Cl. See Tabled for definitions of the listed quantities. 










